NUMERICAL VARIATIONAL METHODS APPLIED TO CYLINDER 
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Abstract. We review and compare different computational variational methods applied to a system of fourth 
order equations that arises as a model of cylinder buckling. We describe both the discretization and implementation, 
in particular how to deal with a 1 dimensional null space. We show that we can construct many different solutions 
from a complex energy surface. We examine numerically convergence in the spatial discretization and in the domain 
size. Finally we give a physical interpretation of some of the solutions found. 

1. Introduction. We describe complementary approaches to finding solutions of systems of 
fourth order elliptic PDEs. The techniques are applied to a problem that arises in the classic 
treatment of an isotropic cylindrical shell under axial compression but are also applicable to a 
wide range of problems such as waves on a suspension bridge [4, 6], the Fucik spectrum of the 
Laplacian [7], or the formation of microstructure [12, 3]. 

The cylindrical shell offers a computationally challenging and physically relevant problem with 
a complex energy surface. We take as our model for the shell the Von Karman-Donnell equations 
which can be rescaled [5] to the form 

A^^ + A^^^ - (^xx - 2 [w, (j)] = 0, (1.1) 
A^0 + ^^^ + [w,w\ = 0, (1.2) 



where the bracket is defined as 



The function w is di scaled inward radial displacement measured from the unbuckled (fundamental) 
state, (j) is the Airy stress function, and A G (0, 2) is a load parameter. The unknowns w and (j) 
are defined on a two-dimensional spatial domain ft = (—a, a) x {—b^b)^ where x G (—a, a) is 
the axial and y G {—b^b) is the tangential coordinate. Since the ^/-domain {—b^b) represents the 
circumference of the cylinder, the following boundary conditions are prescribed: 

w is periodic in and Wx = {^w)x = at x = ±a, (1.4a) 
(j) is periodic in y^ and (j)x = {^(fyx = at x = ±a, (1.4b) 

as shown in Fig. 1.1 (i), (ii). 

1.1. Functional setting. We search for weak solutions (j) of (1.1-1.4) in the space 

X = \^ G H'^{Q) : ipxi^Ci^ ') = is periodic in and J ip = 0^ 

with norm 

||«;||^= / (Aw' + A^l), 
Jn 

where 0i G H'^{^) is the unique solution of 

A^^i = —Wxx^ ^1 satisfies (1.4b), and / ^1 = 0. (1.5) 
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Fig. 1.1. (i) The geometry of the cylinder, (ii) the computational domain and the boundary conditions, (Hi) 
one quarter of the domain and the corresponding boundary conditions. 



This norm is equivalent to the i^^-norm on X, and with the appropriate inner product (•, •)x the 
space X is a Hilbert space. Alternatively, if the load parameter A G (0, 2) is fixed, another norm 

Jn 

can be used. Because of the estimate 

/ wl = - [ WW,, = [ wA^^i = [ AwA^i <l [ Aw^ [ A^l = l\\w\\l, 

JQ JQ JQ JQ ^ JQ ^ JQ ^ 

it is equivalent to and hence also to the i^^-norm on X. The corresponding inner product 

will be denoted (•, •)x,a- 

Equations (1.1-1.2) are related to the stored energy the average axial shortening and 
the total potential given by 

Eiw) ■.= \j (Au;2 + A(/.2) , S{w) ■■=\ f ^1, Fx = E-XS. (1.6) 

Note that the function (f) in (1.6) is determined from w by solving (1.2) with boundary condi- 
tions (1.4b). All the functionals and Fx belong to C^(X), i.e., are continuously Frechet 
differentiable. 

The fact that (1.1) is a reformulation of the stationarity condition = E' — XS' = will be 
important in Sec. 2, and we therefore briefly sketch the argument. It is easy to see that 

S'{w) ■h = - Wxxh. 

JQ 

For E'{w) • h, let w,(j) e X solve (1.2) and h^ip e X solve A'^ip = -hxx - [h, h] - 2[w, h]. Then, 
assuming sufficient regularity on 

E{w + /i) - E{w) = I Aw Ah ^ \ I {Ahf + / A(^AV^ + ^ / ^^^f 

= [ {hA^w-h<P,^-2[w,h]cP) + l [ {Ahf + l [ {A4>f- [ [h,h]<p, 
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where we used integration by parts several times. The last three integrals are for \\h\\^ 

and it can be shown by integration by parts that 

[[w,h]^= [ h[w,^] . (1.7) 

Therefore 

F^{w) -h = E\w) 'h- XS\w) ■h= / /i(^A^^ + A^^^ - (^^^ - 2 0]^ 

1.2. Review of some variational numerical methods. We now describe the variational 
methods used to find numerical approximations of critical points of the total potential Fx- In our 
numerical experiments these methods are accompanied by Newton's method and continuation. 
The advantage of this approach is that it combines the knowledge of global features of the energy 
landscape with local ones of a neighborhood of a critical point. The details related to spatial 
discretization will be discussed in Sec. 2, the Newton-based methods in Sec. 3. 

1.2.1. Steepest descent method (SDM). Let the load parameter A G (0,2) be fixed; 
we work in a discretized version of (X, (•, ■)x,x)- We try to minimize the total potential Fx by 
following its gradient fiow. We solve the initial value problem 

j^w{t) = -VxFx{w{t)) , ^(0) = ^0 , 

with a suitable starting point wq on some interval (0,T]. This problem is then discretized in t. 

In [5] it was shown that = is a local minimizer of Fx. Indeed, if ||i^o|Ixa small, the 
numerical solution w{t) converges to zero as t tends to infinity. If, on the other hand, ||t^^o|lxA 
is large, the numerical solution w{t) stays bounded away from zero. In most of our experiments, 
the numerical algorithm did not converge for t ^ oo in the large norm case. The only exception 
for a relatively small value of A will be mentioned later in Sec. 5.3. Nevertheless, for a sufficiently 
large computational domain ft and a sufficiently large t > we obtain Fx{w{t)) < 0. Such a state 
w{t) is needed for the mountain pass algorithm as explained below. Existence of this state was 
also proved in [5]. 

1.2.2. Mountain-pass algorithm (MPA). The algorithm was first proposed in [2] for a 
second order elliptic problem in ID and extended in [6] to a fourth-order problem in 2D. We give 
a brief description of the algorithm here. 

Let the load A G (0,2) be fixed; we work again in a discretized version of (X, (•, ■)x,x)- We 
denote wi = the local minimum of Fx and take a point W2 such that Fx{w2) < (in practice 
this point is found using the SDM). We take a discretized path {zm}^^Q connecting zq = wi 
with Zp = W2. After finding the point Zm at which Fx is maximal along the path, this point is 
moved a small distance in the direction of the steepest descent —\/xFx{zm)' Thus the path has 
been deformed and the maximum of Fx lowered. This deforming of the path is repeated until the 
maximum along the path cannot be lowered any more: a critical point w-^p has been reached. 
Figure 1.2 illustrates the main idea of the method. 

The mountain-pass algorithm is local in its nature. The numerical solution Wy^p it finds has 
the mountain-pass property in a certain neighborhood only. The choice of the path endpoint W2 
may infiuence to which critical point the algorithm converges. Different choices of W2 are in turn 
achieved by choosing different initial points wq in the SDM. 

1.2.3. Constrained steepest descent method (CSDM). We fix the amount of shorten- 
ing S of the cylinder. This is often considered as what actually occurs in experiments. We work 
now in a discretized version of (A, (•, •)x)- Let C > be a fixed number and define a set of w 
with constant shortening 



^ = {weX : S{w) = C} . 
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(1.8) 



Fig. 1.2. Deforming the path in the main loop of the mountain pass algorithm: point Zm is moved a small 
distance in the direction —'V xFx(zm) CLnd becomes z^"^ . This step is repeated until the mountain pass point w-^p 
is reached. 



Critical points of E under this constraint are critical points of Fx, where A is a Lagrange multiplier. 
The simplest such points are local minima of the stored energy E on We need to follow the 
gradient flow of E on hence we solve the initial value problem 

^w{t) = -P^^t)^E{w{t)) , ^(0) =woe^, 
for t > 0. Pw denotes the orthogonal projection in X on the tangent space of ^ dit w e 

P U = U 7y—\/S(w) . 

\\VS{w)\\l 

The details of the algorithm can be found in [4]. The initial value problem is solved by repeating 
the following two steps: given a point w G ^ find w = w — /SiP^V E{w)^ where At > is small, 
and define it^new = ciD, where the scaling coefficient c is chosen so that ^^;new ^ ^ - The algorithm 
is stopped when \\Py^VE{w)\\^ is smaller than a prescribed tolerance. The corresponding load is 
given by 

{VS{w),VE{w))x 
\\VS{w)\\\ 

1,2 A, Constrained mountain-pass algorithm (CMPA). Let C > and be the set 

of w with constant shortening given in (1.8). We would like to find mountain-pass points of E 
on The method has been described in [4] in detail. We need two local minima wi^W2 of E 
on ^ which can be found using the CSDM. The algorithm is then similar to the MPA. We take 
a discretized path {zm}^=o ^ ^ connecting = vj\ with Zrp = After finding the point 
Zm at which E is maximal along the path, this point is moved a small distance in the tangent 
space to at Zra in the direction of the steepest descent —Pz^VE{^Zm) and than scaled (as in 
the CSDM) to come back to ^ . Thus the path has been deformed on ^ and the maximum of 
E lowered. This deforming of the path is repeated until the maximum along the path cannot be 
lowered any more: a mountain-pass point of E on ^ has been reached. The load A is computed 
as in the CSDM. 

The choice of the end points w\ and W2 will in general infiuence to which critical point the 
algorithm converges. 



1.3. Computational Domains. We consider the problem on the domain Q (Figure 1.1 (ii)) 
both without further restraints and under a symmetry assumption, which reduces the computa- 
tional complexity. In the latter case we assume 



w{x, y) = w{-x, y) = w{x, -y) (1 g) 

(}){x,y) = (l){-x,y) = (l){x,-y) ^ ' 

By looking for solutions w^(j) e X that satisfy (1.9) the domain is effectively reduced to one 
quarter: Qi = (— a,0) x (—6,0) as shown in Figure 1.1 (iii). One needs to solve (1.1-1.2) only on 
Qi with the boundary conditions 

= (A^)^ = (l)u = (A(^)^ = on dQi , (1.10) 

where u denotes the outward normal direction to the boundary. Hence we search for weak solutions 
of (1.1-1.2), (1.10) in the space 



Xi = {il; e H^{ni) : = on 8^1, and / = 



f 



We can then use (1.9) to extend these functions to the whole Q. 

We have performed numerical experiments both with and without the symmetry assumption. 
For the sake of simplicity we will give a detailed description of the numerical methods for the 
second case only where the boundary conditions are the same on all sides ofQi. The first case with 
periodic conditions on two sides of Q is very similar and will be briefly mentioned in Remark 2.1. 

1.4. Solving the biharmonic equation. In order to obtain (j) for a given one has to 
solve (1.2); to compute the norm of one has to solve (1.5). Both problems are of the form 

A'^il; = f inQi, i/j^ = (AiIj)^ = on OQi, [ i/j = Q, (1.11) 

Jn 1 
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where / G L^{Qi) is given. If J^^ / = 0, then (1.11) has a unique weak solution i/j in Xi. It is a 

4 

straightforward calculation to verify that the right-hand sides of equations in (1.2) and (1.5) have 
zero average. 

In the discretization described below the problem (1.11) is treated as a system: 



"^^ ^ in 1^1, = = on Sl^i, / u= I v = 0. (1.12) 



5^1 

4 



/ " 



The system has a unique weak solution (ix, v) G {H^{Qi))'^. Since the domain Qi has no reentrant 
corners. Theorem 1.4.5 of [9] guarantees that v G H'^{Qi) and therefore that the two formulations 
are equivalent. 

2. Finite difference discretization. We discretize the domain ^li by a uniform mesh 
{xm^Vn) ^ with M points in the x-direction and TV points in the ^/-direction: 

= -a + (m - ^)Ax, m G {1 . . . , M}, 

yn = -b^{n- ^)Ay, ne{l...,N}, 

where Ax = a/M, Ay = b/N. We represent the values of some function w on Qi at these points 
by a vector w = {wi)ffj^ ^ where Wi = w{xm^yn) and i = {n — 1)M -\- m. In our notation we will 
not distinguish between as a function and as a corresponding vector. The vector w can also 
be interpreted as a block vector with N blocks, each containing M values of a single row of the 
mesh. We introduce the following conventions for notation: 
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• For two matrices = (a^^O^^^, = (bki)^/^^ we define 
which is an TV X A/" block matrix, each block is an M x M matrix. 

• For two vectors u = {ui)f£^ , v = {vi)f^ we define u Q v = {uiVi)f^ ^ i.e., a product of 
the components. 

To discretize second derivatives we use the standard central finite differences (with Neumann 
boundary conditions [10]). Let Id^ denote the M x M identity matrix and define another M x M 
matrix 



A. 



M 



-1 



The second derivatives —dxx^ ~^yy ^^e biharmonic operator with the appropriate boundary 
conditions are approximated by 

1 



A = 



aM , 



■N 



A = 

yy 



1 



Ay' 



:Id 



'M 



AT 



respectively. 

2.1. Discretization of and the bracket [•,•]. Supposing that we can solve the 

discretized version of (1.2) 

A^2(j) - A^^w + [w, w]2 = 0, (2.1) 
we can also evaluate the energy E and the shortening S: 

E{w) = 2 {w^A^2 w + (/)^ (/)) Ax Ay, S{w) = 2 {w^ A^^w) Ax Ay. (2.2) 

In order to solve (2.1) we need to be able to solve the biharmonic equation and we need to 
choose a discretization of the bracket [•,•]. This bracket appears in the equations in two different 
roles: in equation (1.2) the bracket is part of the mapping w (j), and therefore of the definition 
of the energy E; in equation (1.1), which represents the stationarity condition E' — XS' = 0, the 
bracket appears as a result of differentiating E with respect to w and applying partial integration. 
As a result, we need to use two different forms of discretization for the two cases. 

In both cases the bracket requires a discretization of the mixed derivative dxy One choice is 
to use one-sided finite differences. Define M x M matrices 

■ -1 1 



iM 



1 



-1 1 



iM 



(2.3) 



We choose either left or right-sided differences represented by these matrices, respectively, let Af^ 
denote our choice (cf. Sec. 5.1). The derivatives dx^ dy, and —dxy are approximated by 

(2.4) 



1 

Ax 



y ' 



For the definition of ((> in terms of w (equation (1.2)) we choose 
and the corresponding choice for equation (1.1) is 

[W, <P]l = \Ayy {{A^M © ^} + \Kx {{\yW) © ^} - ^xy {(^xyW^) © 0} 



(2.5) 



(2.6) 



These two definitions are related in the sense given in (1.7): \w ,h}^(\) = In^ \w , iox all h. 

With these definitions the partial derivatives of discretized E and S with respect to the 
components of w are given by 



E'(w) = A^2W + A^ 



2[^,( 



Ji 5 



S'(w) 



^xx^' 



(2.7) 
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2.2. Solving the discretized biharmonic equation. Matrix A^2 is symmetric and has 
a one-dimensional nullspace: A^2l = 0, where 1 and are vectors with MA^-components which 



are all one and all zero, respectively. The same is true for 
would like to solve 



and Ayy. For a given vector / we 



A^2tlj = f, 



0. 



(2.8) 



A unique solution exists if and only if / has zero average, i.e.. I f = 0. So we must verify that 
the discretized versions of the right-hand sides in (1.2), (1.5) satisfy this condition. Let w G R^^, 
then 



1 A^x^ 



0, 



l^[W, = {A^^wf{AyyW) - {A^yWf{A^yW) = A ^ ^ A y y W " A^ y A ^ y W = , (2.9) 



where the last equality holds because A^A^ = A^^ and Ay Ay = Ayy^ and because the x- and 



^/-matrices commute. We have, in fact, shown that the integration by parts formula from the 
continuous case holds for our choice of spatial discretization. This is not true for an arbitrary 
discretization but is key for a successful scheme. 

The inverse of matrix A^2 on the subspace of vectors with zero average, denoted with a slight 
abuse of notation by A^2, can be found, for example, using the fast cosine transform described 
below in Sec. 2.4. 

2.3. Computing the gradient. The variational methods of this paper are based on a steep- 
est descent flow and modifications of this algorithm. The direction of the steepest descent of E at 
a point w G X is opposite to the gradient of E dit w. The gradient is the Riesz representative of 
the Frechet derivative and hence we need to find a vector u e X, such that E'{w) - v = {u, v) for 
all V e X. The inner product is either (•, •)x or (•, -)x,x and hence the gradient depends on the 
choice of the inner product. We use the notation u = \/E{w) for the gradient in (X, (•, •)x) and 
u = \/xE{w) for the gradient in (X, (•, ■)x,x)- To find the discretized version of the gradient, we 
first need to discretize the inner product. 

Let u^v ^ R^^, l^u = l^v = 0. The inner product is evaluated in the following way: 



{u, v)x,x = 4 (u^A^2V + (l)f^A^2(l)l - Xu^A^^v^ Ax /\y 
= a(u^ {a ^2 



A A 

^xx^A'^^xx 



are solutions of the discretized version of (1.5) with w replaced by u and respec- 



where (f)\ 

tively, and we assume that we work on Qi. For w G 
E'{w) given in (2.7) is computed as 



the Riesz representative of 



A A~^ A 

^xx^A"^ XX 



E\w) . 



(2.10) 



As in the case of we abused notation here since the inverse only makes sense on a subspace of 
vectors with zero average. It can be easily verified that 1^ E'{w) = 0. The numerical evaluation 
of V \S and of the (•, •)x-gi'adients is similar. 

2.4. Fourier coordinates. In Fourier coordinates most of the finite difference operators 
become diagonal matrices. This increases the efficiency of the numerical algorithm and makes it 
possible to easily find the inverse of matrices like A^2. See for example [1]. 

On a uniform grid, it is a standard procedure to apply some form of the fast Fourier transform 
to diagonalize finite difference matrices like (see, e.g., [11]). Due to the Neumann boundary 
conditions (1.10) we need to employ the fast cosine transform. We define M x M matrices 



/2M 



2 cos 



(i-l)(2j-l)^ 
2M 



M 



/2M 



2 cos 



(2i-l)(j-l)7r 
2M 



M,M 



^=l,j=2 
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which have the fohowing properties: 



^2 '-^b 



A 



M 



where A 



M 



AM 
^2 



We further define matrices 

C — 

which diagonahze A^^, Ayy^ and A ^2'. 

where the diagonal matrices are given by 

A_ = ^A^®Jd^, A, 
For a vector w G 



2 cos Hence they are inverses of each other and diagonahze 



f 5 



b •> 



C^AyyC^ 



C^A^2C^ — A^2, 



yy 



i A 



AT 



A )2. 



(2.11) 
(2.12) 



we introduce its Fourier coordinates w by 

w = C^w w = C^^w . 

We note that l^w = if and only if the first component of w is zero. 

Most of the computations involved in the variational methods described in Sec. 1.2 can be done 
in the Fourier coordinates. The only time one needs to go back to the original coordinates is when 
evaluating the brackets (2.5) and (2.6), because they are nonlinear and involve the discretized 
mixed derivative operator A^y. 

2.5. Alternative discretization of —dxy The fast Fourier transform provides us with 
another discretization of the mixed derivative which is not biased to the left or right. In an 
analogy to (2.12) and (2.11) we define 



xy Ax Ay 

where S is the fast sine transform matrix 



A^0a/A^, 



'■xy 



^■^xy^i 



S = S 



M , 



J AT 



s 



M 



/2M 







(2 sin 



(2i-l)(i-l) 
2M 



M,M 
,J=2 



Property (2.9) also holds with A^y replaced by A^y. 

Remark 2.1. When discretizing the problem on the full domain Q with boundary conditions (1.4), 
we need to use different matrices in the x and ^/-directions. In the x-direction we use the matrices 
described above, in the ^/-direction to discretize the second derivatives, for example, we use 



-1 



In this direction the fast Fourier transform is used instead of the fast cosine/sine transform. 

3. Newton's method. We use Newton's method in two different ways. The first is to 
improve the numerical approximations obtained by the variational numerical methods. Since 
these are sometimes slow to converge, it is often faster to stop such an algorithm early and use its 
result as an initial guess for Newton's method. The second use for Newton's method is as part of 
a numerical continuation algorithm (see Sec. 3.3). 



3.1. Newton's method for given load parameter A. This method can be used to improve 
solutions obtained by the MPA. Let A G (0, 2) be given. We are solving 





' ^1 " 




A^2W - XA^^w + A^^cf) - 2[w, 0]i 




" " 




. ^2 _ 




-A^2(t) + A^^w - [w, w]2 








(3.1) 



for w and (j) with zero average using Newton's method. The matrix we need to invert is 









dw 








3^2 




- dw 


del) . 





where 



d 1 

5 r n 1 5 r 



2Bi 



2Bi 



- 2-B2 



(3.2) 



T 



^yy(diag( 



(diag(/))A 



-A^^(diag A^^^) - A^^(diag A^^^) . 



From the properties of A^^, A , A^ , A^2 it follows that matrix W{w^(j)) is symmetric and sin- 



gular, its nullspace is spanned by 



" 1 







1 



To describe how to find the inverse of ^'{w^ (j)) on a subspace orthogonal to its nullspace, we 
introduce a new notation for the four blocks of ^'{w^ (j)) from (3.2): 

For given vectors li, r] with zero average we need to find vectors ( with zero average such that 

(3.3) 





^12 




V 




u 


^12 


^22 




. c . 







Let the tilde denote the block of the first MN — 1 rows and columns of a matrix and MN — 1 

^11 ^12 



components of a vector. The matrix 



^12 ^22 



is symmetric, nonsingular, and sparse. It can 



be inverted by a linear sparse solver. System (3.3) is then solved in the following steps: 











-1 




r 






G12 




u 














. P _ 




G12 


G22 







V 



r -\- si 
s 

(J 



3.2. Newton's method for given S. This method can be used to improve solutions ob- 
tained by the CSDM and the CMPA. Let C > be given. We are looking for numerical solutions 
of (1.1-1.2) in the set ^ defined by (1.8). Hence we are solving 



(3.4) 



for w and (/) with zero average and A using Newton's method. The approach is very similar to that 
described in the previous section, the resulting matrix is symmetric, has just one more row and 
column than the matrix in (3.2). 

9 



A^2W - XA^^w + A^^(/) - 2[w, 




' " 


- A^2(/) + A^^w — [w, w]2 







-^w^A^^w^C/{4AxAy) 








3.3. Continuation. To follow branches of solutions (X^w) of (3.1) we adopt a continuation 
method described in [8]. We introduce a parameter s G M by adding a constraint — pseudo- 
arclength normalization (in the (A, j^)-plane). For a given value of s we are solving 





" ^1 " 




A^2W - XA^^w + - 2[w, (^]i 




' " 








-A^20 + A^^w - [w, w]2 









_ % _ 




_ 9{wo,w - wq)x + (1 - 0)Xo{X - Ao) - {s - sq) 








(3.5) 



for (j) with zero average, and the load A, where the value of ^ G (0, 1) is fixed (e.g., = 
^). We assume that we are given a value sq, an initial point (Ao,i^o) on the branch, and an 
approximate direction {Xq^wq) of the branch at this point (an approximation of the derivative 

System (3.5) is solved for a discrete set of values of s in some interval (5o,5i) by Newton's 
method. Then, a new initial point on the branch is defined by setting wq = w{si)^ Aq = A(5i), 
So = 5i, a new direction (Xq^wq) at this point is computed and the process is repeated. The 
matrix we need to invert in Newton's method is 



P(^,(/),A) 



^'(^,(/),A) g 
d 







AOA^^^^WQ Ax Ay 




A A~IA 

^xx^A^ XX' 



d=il-0)Xo, 
(3.6) 



where = A^2 

Solving a linear system with this matrix amounts to solving system (3.3) for two right-hand 
sides. For a given u G M?^^ with [1-^ l^]u = and a given 77 G M we want to solve 



\w,(t),X) g 




V 




u 


d 




. c . 







(3.7) 



for V with [1 1 ]v = and (. System (3.7) is solved in the following steps: 



solve : ^\w^ (/), A)'^! = g , 



d — h^vi ' 



V = V2 — (Vi . 



Remark 3.1. Note that in this implementation we simply follow a solution of the equation, there 
is no guarantee that this remains a local minimum, a MP-solution or constrained MP-solution 
(cf. Fig. 5.6). 



Remark 3.2. Newton's method and continuation have been implemented only using a one-sided 
discretization of the mixed derivative dxy and only on the domain assuming symmetry (1.9). 
The alternative discretization of dxy described in Sec. 2.5 uses the fast cosine/sine transform. The 
resulting matrix A^y is not sparse and therefore we would obtain a dense block G12 in system (3.3) 
which would prevent us from using a sparse solver. 

On the full domain we assume periodicity of w and (j) in the ^/-direction. Hence for a 
discretization with a small step Ay the matrix we invert when solving (3.3) would become close to 
singular. The shift in the y direction is prevented by assuming the symmetry w{x,y) = w{x^ —y), 
(t){x,y) = (t){x,-y). 

4. Numerical solutions. We fix the size of the domain and the size of the space step for the 
following numerical computations: a = 6 = 100, Ax — Ay = 0.5. We obtain solutions using the 
variational techniques SDM, MPA, CSDM, and CMPA. Table 4.1 provides a summary of which 
discretization was used in which algorithm. 
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mixed derivative dxy 


computational domain 


one-sided 


Fourier 


full 


1/4 


variational methods 


^/ 


^/ 


^/ 


^/ 


Newton/continuation 











Table 4.1 

Summary of which spatial discretization was used in the different numerical methods. 



4.1. A mountain-pass solution on the full domain ft. The first numerical experiments 
are done on the full domain 1], without symmetry restrictions, and with the unbiased (Fourier) 
discretization of dxy (Sec. 2.5). For a fixed load A = 1.4 we computed a mountain-pass solution 
using the MPA (Sec. 1.2.2). As end points were taken wi = and a second point W2 obtained 
by the SDM (here the initial point for the SDM was chosen to have a single peak centered at 
X = y = 0). The graph of the solution Wy^p is shown in Fig. 4.1 (left). The figure on the right 
shows w^p rendered on a cylinder and we see it has the form of a single dimple. The value of 
shortening for this solution is S{wy^p) = 14.93529. 

Alternatively, if we apply the CSDM with S = 14.93529 and use a function with a single peak 
in the center of the domain as the initial condition wq we also obtain the same solution w-^p, this 
time as a local minimizer of E under constrained S. 

We remark that although the MPA and the CSDM have a local character we have not found any 
numerical mountain-pass solution with the total potential Fx smaller than Fx{wy^p) for A = 1.4. 
Similarly, using the CSDM we have not found any solution with energy E smaller than E{wy^p) 
under the constraint S = 14.93529. We briefly discuss the physical relevance of this solution 
in Sec. 6. 




Fig. 4.1. Mountain-pass solution for X = 1.4 found using the MPA on the full domain Q with dxy discretized 
using the fast Fourier transform. 

4.2. Solutions under symmetry restrictions. The solution Wy^p of Fig. 4.1 satisfies the 
symmetry property (1.9). In the computations described below we enforced this symmetry and 
worked on the quarter domain l^i, thus reducing the complexity of the problem. In order to 
improve the variational methods by combining them with Newton's method we also discretized 
the mixed derivative dxy using left-sided finite differences. The influence of this choice on the 
numerical solution is described in Sec. 5.1. 

4.2.1. Constrained steepest descent method. We first fixed /S = 40 and used the CSDM 
to obtain constrained local minimizers of E described in Table 4.2. They are ordered according 
to the increasing value of stored energy E. Their graphs and renderings on a cylinder are shown 
in Fig. 4.2. Solution 1.1 is similar to the single dimple solution Wy^p described above and according 
to Table 4.2 it has, indeed, the smallest value of E. 
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CSDM 


A 


S 


E 


Fx 


same shape as MPA 


1.1 


1.108121 


40 


56.85636 


12.53151 


2.1 


1.2 


1.299143 


40 


62.76150 


10.79577 


2.2 


1.3 


1.316146 


40 


63.21646 


10.57063 


2.3 


1.4 


1.311687 


40 


63.64083 


11.17334 


2.4 


1.5 


1.309586 


40 


63.70623 


11.32278 


2.5 


1.6 


1.328997 


40 


64.00875 


10.84889 


2.6 


1.7 


1.344898 


40 


64.52244 


10.72651 


2.7 



Table 4.2 

Numerical solutions obtained by the CSDM on Qi with dxy discretized using left-sided finite differences. 

4 

Graphs are shown in Fig. 4-^- 
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Fig. 4.2. Numerical solutions found using the CSDM with 



axial end shortening S = 40. More details are in Table J^..2. 



4.2.2. Mountain-pass algorithm. We then used the MPA for fixed A = 1.4 and various 
choices of W2 to obtain the local mountain-pass points of Fx described in Table 4.3. They are 
ordered according to the increasing value of the total potential Fx. The shape of their graph is 
very similar to that of the CSDM solutions discussed above and depicted in Fig. 4.2 and we do 
not show their graphs here. Solution 2.1 is again the single dimple solution and the table shows 
that it has the smallest value of Fx. 



MPA 


A 


S 


E 


Fx 


same shape as CSDM 


2.1 


1.4 


17.73822 


29.42997 


4.596460 


1.1 


2.2 


1.4 


29.85121 


49.08882 


7.297132 


1.2 


2.3 


1.4 


31.28849 


51.39952 


7.595635 


1.3 


2.4 


1.4 


31.41723 


52.01893 


8.034809 


1.4 


2.5 


1.4 


31.22992 


51.84074 


8.118852 


1.5 


2.6 


1.4 


32.77491 


54.15818 


8.273314 


1.6 


2.7 


1.4 


34.19888 


56.56472 


8.686284 


1.7 



Table 4.3 

Numerical solutions obtained by the MPA onfli with dxy discretized using left- sided finite differences. 



4.2.3. Constrained mountain-pass algorithm. We then fixed 6* = 40 and applied the 
CMPA to obtain constrained local mountain passes of E described in Table 4.4. They are again 
ordered according to the increasing value of stored energy E. Their graphs and renderings on a 
cylinder are shown in Figs. 4.3 and 4.4. As endpoints W2 of the path in the CMPA we used 
the constrained local minimizers 1.1-1.7. 

There are 21 possible pairs {wi^W2) to be used but only 19 solutions in Table 4.4. The 
algorithm did not converge for the following three pairs: (1.1, 1.3), (1.5, 1.6), and (1.4, 1.6), 
most likely due to the complicated nature of the energy landscape between these endpoints. On 
the other hand, two choices of pairs denoted by * and f in the Table yielded two solutions each. 
When the path is deformed it sometimes comes close to another critical point of E which is not a 
constrained mountain pass. In that case the algorithm slows down and one can apply Newton's 
method to such a point. It is a matter of luck whether Newton's method converges. The CMPA 
then runs further and might converge to another point, this time a constrained mountain-pass 
point. And finally, two choices of {wi^W2) yielded the same solution 3.3. 
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CMPA 


A 


S 


E 


Fx 


end points 


3.1 


1.310815 


40 


63.98996 


11.55737 


1.2, 1.4 


3.2 


1.332112 


40 


64.38609 


11.10161 


1.3, 1.6 


3.3 


1.447626 


40 


66.49032 


8.585294 


1.2, 1.3 or 1.3, 1.4 


3.4 


1.440841 


40 


66.72079 


9.087129 


1.1, 1.4 


3.5 


1.447594 


40 


66.97057 


9.066810 


1.1, 1.6 


3.6 


1.484790 


40 


68.20637 


8.814758 


1.3, 1.5 


3.7 


1.477769 


40 


68.23274 


9.121955 


1.1, 1.5* 


3.8 


1.482261 


40 


68.41086 


9.120428 


1.1, 1.7 


3.9 


1.413917 


40 


68.56697 


12.01028 


1.1, 1.2 


3.10 


1.520975 


40 


68.83818 


7.999162 


1.2, 1.5 


3.11 


1.475705 


40 


69.00087 


9.972652 


1.1, 1.5* 


3.12 


1.532000 


40 


69.27379 


7.993781 


1.4, 1.5"^ 


3.13 


1.527108 


40 


69.35834 


8.274019 


1.2, 1.6 


3.14 


1.551762 


40 


69.47838 


7.407904 


1.4, 1.5"^ 


3.15 


1.547955 


40 


69.68292 


7.764712 


1.6, 1.7 


3.16 


1.539785 


40 


69.78487 


8.193480 


1.5, 1.7 


3.17 


1.546480 


40 


69.85900 


7.999795 


1.3, 1.7 


3.18 


1.549780 


40 


70.16253 


8.171339 


1.2, 1.7 


3.19 


1.561117 


40 


70.74117 


8.296474 


1.4, 1.7 



Table 4.4 

Numerical solutions obtained by the CMPA /Newton onfli with dxy discretized using left- sided finite differ- 

4 

ences. Graphs are shown in Figs. 4-3, 4-4- 
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Fig. 4.3. Numerical solutions found using the CMPA/Newton with S = 40. More details are in Table 4-4- 




Fig. 4.4. Numerical solutions found using the CMPA/Newton with S = 40. More details are in Table 4-4- 



5. Remarks on the numerics. 



5.1. Bias in the discretization of dxy In this section we examine the influence of the 
discretization of the mixed derivative dxy on the numerical solution. We recall that the mixed 
derivative dxy can be discretized using left /right-sided finite differences (2.3), (2.4) or using the 
fast Fourier transform (Section 2.5). For comparison we use the single-dimple solution on 1] = 
(-100, 100)2 at load A = 1.4 obtained by the MPA. 

Let Ax = = 0.5. Table 5.1 gives a list of numerical experiments together with the values of 
shortening and energy. Figure 5.1 shows a profile of the numerical solutions in the circumferential 
direction at x = 0. 



domain 


discretization of dxy 


A 


S 


E 


Fx 


Figure 5.1 


1] or 1^1 

4 


Fourier 


1.4 


14.93529 


24.71825 


3.808850 






left/right-sided 


1.4 


14.93617 


24.70828 


3.797636 






left-sided 


1.4 


17.73822 


29.42997 


4.596460 




4 


right-sided 


1.4 


12.81205 


21.16342 


3.226549 





Table 5.1 

Single-dimple numerical solution obtained by the MPA with and without the symmetry assumption (1.9) and 
with various kinds of discretization of dxy- 




On the full domain with no assumption on symmetry of solutions the discretization of dxy 
using the left/right-sided finite differences (2.3) provides a numerical solution that is slightly asym- 
metric (Fig. 5.1, thin solid line). The Fourier transform provides a symmetric solution (Fig. 5.1, 
thick solid line). The same numerical solution can be obtained on ^li under the symmetry as- 
sumption (1.9) with dxy discretized using the fast cosine/sine transform. 

On Vti the symmetry of numerical solutions is guaranteed by assumption (1.9). The use of 
left/right-sided discretization of dxy does, however, have an influence on the shape of the numerical 
solution, as Fig. 5.1 shows (the dotted and the dashed line). 

5.2. Convergence. We now turn to the influence of the size of the space step Ax, A^ on 
the numerical solution. We run the MPA on Vti under the symmetry assumption (1.9) with dxy 
discretized using (a) the fast cosine/sine transform, (b) left-sided finite differences, (c) right-sided 
finite differences. We consider Ax = A?/ = 0.5, 0.4, 0.3, 0.2, 0.1, i.e., we take 200, 250, 333, 500, and 
1000 points in both axis directions, respectively. Figure 5.2 illustrates convergence as Ax, A?/ 
of the numerical solutions obtained by various types of discretization of dxy 
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0.2 



0.3 
A x=A y 



0.4 



0.5 
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16 



13 
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S{wcs) 
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0.2 0.3 0.4 

Ax=Ay 



Fig. 5.2. Influence of the size of the space step Ax, Ay on the numerical solution wmp obtained by the MPA 
for three different kinds of discretization of dxy Let Wj^, denote the numerical solutions obtained using the 
left and right-sided discretization of dxy, respectively, w^g using the fast cosine/sine transform. Left: comparison 
of the solutions in the maximum norm; right: the value of shortening S. 



5.3. Dependence on the size of the domain. As observed in [5], the locaHzed nature 
of the solutions suggests that they should be independent of domain size, in the sense that for a 
sequence of domains of increasing size the solutions converge (for instance uniformly on compact 
subsets). Such a convergence would also imply convergence of the associated energy levels. Sim- 
ilarly, we would expect that the aspect ratio of the domain is of little importance in the limit of 
large domains. 

We tested these hypotheses by computing the single-dimple solution on domains of different 
sizes and aspect ratios. In all the computations the space step Ax = A?/ = 0.5 is fixed. In order to 
use the continuation method of Sec. 3.3, we discretized dxy using the left-sided finite differences. 
We also assumed symmetry of solutions given by (1.9) and worked on Vti. We recall the notation 
of computational domains, Vt = (—a, a) x (—6,6), l^i = (— a,0) x (—6,0). 

Figure 5.3 shows the results for load A = 1.4. First we notice that the central dimple has 
almost the same shape in all the shown cases. But there seems to be a difference in the slope of 
the "flat" part leading to this dimple. On domains with small a (short cylinder) the derivative 
in the circumferential ^/-direction in this part is larger than on domains with larger a (longer 
cylinder). The circumferential length h seems to be less important for the shape of the solution: 
for example, the cases (200,50) and (200,100) look like restrictions of the case (200,200) to smaller 
domains. 

We take a closer look at domains of sizes (a, 6) = (100,100), (100,200), (200,100), and 
(200,200) and the corresponding solutions it;ioo,ioo, '^100,200, '^200,100, and 1^^200,200 shown in the 
figure. We compare the first three with the last one, respectively. It does not make sense to 
compare the values of w itself since the energy functional Fa depends on derivatives of w only. 
We choose to compare Wxx and Wyy. Table 5.2 gives the infinity norm of the relative differ- 
ences. Fig. 5.4 shows graphs of the difference ^^;loo,lOO — '^200,200 and of the second derivatives 
(^200,200 - ^ioo,ioo)xx, (^200,200 - ^100,100)2/2/ on the subdomain (-100,0)^. 

We conclude that solutions on different domains compare well; the maximal difference in 
the second derivatives of w is three orders of magnitude smaller than the supremum norm of 
the same derivative. We also observe that varying the length parameter a while keeping the 
circumference parameter h fixed causes larger changes in the numerical solution than varying the 
cylinder circumference while keeping the length fixed. 

Another way of studying the influence of the domain size on the numerical solution is com- 
paring solution branches obtained by continuation as described in Sec. 3.3. We start with the 
mountain-pass solution for A = 1.4 shown in Fig. 5.3 and continue it for both A > 1.4 and A < 1.4. 
The results are presented in Fig. 5.5. We observe that the branches corresponding to the con- 
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(a, b) = (50, 50) (a, b) = (50, 100) (a, 6) = (50, 200) 




Fig. 5.3. The single-dimple mountain-pass solution with A = 1.4 computed under the assumption of symme- 
try (1.9) with left-sided discretization of dxy for various domain sizes. 



'W^100,100 — ^200,200 ("^^200, 200 — '"^100, 100 )ccx ('"^200,200 — '"^100, 100 ) yy 




Fig. 5.4. Comparison of solutions i(;ioo,ioO; 'w;200,200 from Fig. 5.3 obtained on square domains with a = b = 
100 and a = b = 200^ respectively, and their second derivatives. For a reference, we note that \\{w200,2Qq)xx\\^ = 
1.064522 and ||(^200,20o)yy ||^ = 0.8242491. 





\\[^W-W200,200)xx\\^ 


||(w-W200,200)yy||^ 


(W200,200)xx \\^ 


||(w200,200)yy 11^ 


^ = ^100,100 


2.835 • 10-3 


5.313-10-3 


^ = ^100,200 


1.943 • 10-3 


4.917-10-3 


^ = ^200,100 


1.827-10-'^ 


9.638 - 10-"^ 



Table 5.2 

Comparison of the second derivatives of solutions from Fig. 5.3 computed on domains with (a, b) = (100, 100), 
(100,200), (200,100), and (200,200). 
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sidered domains do not differ much for the range of A between approximately 0.71 and 2. Below 
A ~ 0.71 the size of the domain, particularly the length of the cylinder described by a, has a strong 
influence. The graph on the right shows that the larger (longer) the domain Q the smaller the 
value of A at which the norm \\w\\^ starts to rapidly increase for decreasing A. The graph on the 
left shows the energy Fx{w) along a solution branch. The data shown here correspond to the ones 
in the graph on the right marked by a solid line. The dashed line in the right graph shows also 
some data after the first limit point is passed. 




0.5 1 1.5 ,2 Q.5 D.B 0.7 



Fig. 5.5. Continuation of the single-dimple solution found as a numerical mountain pass for A = 1.4 on 
domains of various sizes for a range of values X. Left: Fx(w) as a function of X, right: \\w\\^ as a function of X. 

Figure 5.6 shows how the graph of w{x^y) changes as a solution branch is followed. Here we 
chose a square domain with a = b = 100 and plotted the solution for four values of A (note that 
Figs. 5.6(c), 5.3(100,100), and 5.1(dotted line) show the same numerical solution). We observe 
that with decreasing A the height of the central dimple increases, the dimple becomes wider, and 
the ripples (present at A close to 2) disappear. In Fig. 5.6(a) we observe that new dimples are 
being formed next to the central dimple. 



(a) A ~ 0.593 (b) A = 0.61 (mountain pass) (c) A = 1.4 (mountain pass) 




Fig. 5.6. A detailed look at the continuation of the single- dimple solution on the domain with a = b = 100. 

It should be remarked that although we started the continuation at A = 1.4 at a mountain-pass 
point, not all the points along a continuation branch are mountain passes. Since it is not feasible 
to use the MPA to verify this for each point, we chose just a few. Still on the example of the 
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domain with a = b = 100 in Fig. 5.6, the circles on the continuation branch mark those points that 
have also been found by the MPA (for A = 0.61, 0.65, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 1.95). As described 
in Sec. 1.2.2, in order to start the MP- algorithm a point W2 is needed such that Fx{w2) < 0. The 
analysis in [5] shows that for a given A such a point exists provided the domain Q is large enough 
and in practice it is found by the SDM of Sec. 1.2.1. This was, indeed, the case for all the chosen 
values of A except for A = 0.61. In this case, starting from some wq with a large norm, the SDM 
provides a trajectory w{t) such that Fx{w{t)) > for all t > 0. In fact, the steepest descent 
method converges to a local minimizer w-^ with Fx{wy^) ~ 76.1. This is hence no mountain pass 
but, nevertheless, lies on the same continuation branch and is marked by a triangle in the figure. 
Despite Fx{wy^) > we can still try to run the MPA with W2 = w-^. It converges and yields w-^p 
with Fx{wy^p) ^94.8 (marked by a circle at A = 0.61 and shown in graph (b)). 

The comparison of solutions computed on different domains and their respective energies 
suggests that for each A we are indeed dealing with a single, localized function defined on R^, 
of which our computed solutions are finite-domain adaptations. Based on this suggestion and 
the above discussion of the mountain-pass solutions we could, for example, conclude that the 
mountain-pass energy 

V{X,n) :=mf{Fx{wMp{\n,W2)) :Fa(^2) <0} 

is a finite-domain approximation of a function V^(A), whose graph almost coincides with that of 
V{X,n) for A not too small (cf. Fig. 5.5 left). 

6. Discussion. 

6.1. Variational numerical methods. We have seen that given a complex energy surface 
many solutions may be found using these variational techniques. For example, for a fixed end 
shortening of S = 40, Fig. 4.2, Fig. 4.3 and Fig. 4.4 are all solutions. Which of these solutions is 
of greatest relevance depends on the question that is being asked. 

In the context of the cylinder (and similar structural applications) the mountain-pass solution 
from the unbuckled state {wi = 0) with minimal energy is of physical interest. Often the experi- 
mental buckling load may be at 20-30% of the linear prediction from a bifucation analysis (in our 
scaling this corresponds to A = 2). This uncertainty in the buckling load is a drawback for design 
and so "knockdown" factors have been introduced, based on experimental data. It was argued 
in [5] that the energy of the mountain-pass solution Wy^p in fact provides a lower bound on the 
energy required to buckle the cylinder and so these solutions provide bounds on the (observed) 
buckling load of the cylinder. 

This example illustrates an important aspect of the (constrained) mountain-pass algorithm: 
its explicit non-locality. The algorithm produces a saddle point which has an additional property: 
it is the separating point (and level) between the basins of attraction of the end points wi and W2. 

Another technique to investigate a complex energy surface is to perform a simulated annealing 
computation, essentially to solve the SDM (or the CSDM) problem with additive stochastic forcing. 
The aim in these techniques is often to find a global minimizer (if it exists) where there are a large 
number of local minimizers. Here by either the MPA or the CMPA we find the solution between 
two such minima and so get an estimate on the surplus energy needed to change between local 
minima. 

6.2. Numerical issues. The numerical issues that we encountered are of two types. First 
there are the requirements that are related to the specific problem of the Von Karman-Donnell 
equations, such as the discretization of the mixed derivative and the bracket, and the fact that 
the solutions are symmetric and highly localized. 

For other difficulties it is less clear. For smaller values of A each of the variational methods 
converged remarkably slowly. Newton's method provides a way of improving the convergence, but 
the question is relevant whether this slow convergence is typical for a whole class of variational 
problems. It would be interesting to connect the rate of convergence of, for instance, the SDM to 
certain easily measurable features of the energy landscape. 
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